*Install packages sppack, shp2dta, st0243_1, coefplot, zscore, estout, outreg2
cd "/Users/dima/Dropbox/Slavic Review 1917/FGK Slavic Review Replication"


**************************
*****1917 Unrest DATA*****
**************************

**
**Data from Provisional Government
**

cd Unrest_1917

*Append data for March - September into one dataset
use "Peasant_1917_PG_September.dta", clear
rename gub_map_europe region
gen period = "_September17"
save temp, replace

local month "August July June May April March"
foreach x of local month {
use "Peasant_1917_PG_`x'.dta", clear
rename gub_map_europe region
gen period = "_`x'17"
append using temp
save temp, replace
}

*Prepare disturbance variables to be used in analysis

*Order the dataset
order region period seizure_estate* seizure_land* seizure_harvest* ///
seizure_meadow* seizure_forest* seizure_tools* compulsory* *forest ///
remove* destruct* arson* taxation* other* total rural_population_1916

*Sum up "organized" and "unorganized" disturbances by type if they are 
*given separately
replace seizure_estate = seizure_estate_org + seizure_estate_unorg ///
if seizure_estate == .
replace compulsory_land_rent = compulsory_land_rent_org + ///
compulsory_land_rent_unorg if compulsory_land_rent == .
replace prohib_cut_forest = prohib_cut_forest_org + ///
prohib_cut_forest_unorg if prohib_cut_forest == .
replace remove_work = remove_work_org + remove_work_unorg if remove_work == .
replace other = other_org + other_unorg if other == .

*Other seizures excluding land seizure
egen seizure_rest = rowtotal(seizure_meadow_org seizure_meadow_unorg ///
seizure_forest_org seizure_forest_unorg seizure_tools_org ///
seizure_tools_unorg seizure_harvest_org seizure_harvest_unorg)

*Disruption of activities
egen prevent = rowtotal(prohib_cut_forest remove_work)

*Drop redundant variables
drop prohib_cut_forest remove_work

*Estate destruction
gen destruct = destruct_estate_org + destruct_estate_unorg

*Drop redundant variables
drop *org 

*Recode month variables
gen date = 317 if period == "_March17"
replace date = 417 if period == "_April17"
replace date = 517 if period == "_May17"
replace date = 617 if period == "_June17"
replace date = 717 if period == "_July17"
replace date = 817 if period == "_August17"
replace date = 917 if period == "_September17"

*Drop the previous month variables
drop period

*Reshape the dataset into wide format: 
reshape wide seizure* compulsory_land_rent prevent cut* arson*  destruct ///
other total rural_population_1916, i(region) j(date)

*Rural population as of 1/1/16
rename rural_population_1916317 ruralpopn_1916
drop rural*17
replace ruralpopn_1916 = ruralpopn_1916/100000

*Order the dataset
order region ruralpopn_1916 seizure* compulsory_land_rent* prevent* cut* ///
arson* destruct* other* total*

*March-September total number of disturbances
egen total_1917 = rowtotal(total*)

*April-September totals of disturbances by category
egen seizure_rest_as1917 = rowtotal(seizure_rest*)
gen seizure_estate_as1917 = seizure_estate417+seizure_estate517+ ///
seizure_estate617+seizure_estate717+seizure_estate817+seizure_estate917
gen compulsory_land_rent_as1917 = compulsory_land_rent417+ ///
compulsory_land_rent517+compulsory_land_rent617+compulsory_land_rent717+ ///
compulsory_land_rent817+compulsory_land_rent917
gen prevent_as1917 = prevent417+prevent517+prevent617+prevent717 + ///
prevent817+prevent917
gen other_as1917 = other417+other517+other617+other717+other817+other917

*Disturbances per capita
gen PG_seizure_estate_as1917_pc = seizure_estate_as1917 / ruralpopn_1916
gen PG_seizure_rest_as1917_pc = seizure_rest_as1917 / ruralpopn_1916
gen PG_comp_land_rent_as1917_pc = compulsory_land_rent_as1917 / ruralpopn_1916
gen PG_prevent_as1917_pc = prevent_as1917 / ruralpopn_1916
gen PG_other_as1917_pc = other_as1917 / ruralpopn_1916
gen PG_total_1917_pc = total_1917 / ruralpopn_1916

*Keep disturbance variables and population
keep region ruralpopn PG*
sort region

*Save in file
save ../work, replace

**
**Data from Malyavskii (1981)
**

use "Peasant_1917_Malyavskii.dta", clear
*Note: For Astrakhan' and Saratov only joint total available in Malyavskii's 
*data, so divided the total proportionally to the number of 
*disturbances in the PG data.

*Drop the rural population variable 
drop rural

*Merge with the PG dataset
merge 1:1 region using ../work
drop _m

*Total March - October disturbances per capita
gen MAL_totaltoOct_1917_pc = total_1917 / ruralpopn_1916
gen MAL_totaltoSept_1917_pc = (total_1917 - october) / ruralpopn_1916

*Don't calculate monthly disturbances: for many regions, 
*only have sum across a few months.
*Drop old variables with total and monthly disturbances 
drop total_1917 march-october

*Save in file
save ../work, replace


*************************
*****OTHER VARIABLES*****
*************************

cd ../Other_Variables

*Add 1858 population data
insheet using Bushen_rural.csv, clear
merge 1:1 region using ../work
drop _m
save ../work, replace

*Add data on the number serfs and calculate it per capita
insheet using Bushen.csv, clear
keep name field* house*
mvencode field* house*, mv(.=0) override
gen serfs = fieldserfs_m + fieldserfs_f + houseserfs_m + houseserfs_f
drop field* house*

*Rename "name" variable
rename name region

**Harmonize region names
replace region = "Arkhangel'sk" if region == "Arkhangelsk"
replace region = "Astrakhan'" if region == "Astrakhan"
replace region = "Don" if region == "Don Voisko"
replace region = "Estland" if region == "Estonia"
replace region = "Kazan'" if region == "Kazan"
replace region = "Grodno" if region == "Hrodna"
replace region = "Khar'kov" if region == "Kharkov"
replace region = "Kurland" if region == "Kurliandia"
replace region = "Lifland" if region == "Livonia"
replace region = "Nizhegorod" if region == "Nizhni Novgorod"
replace region = "Perm'" if region == "Perm"
replace region = "Podolsk" if region == "Podolia"
replace region = "Ryazan'" if region == "Ryazan"
replace region = "St.Petersburg" if region == "Petersburg"
replace region = "Tavrida" if region == "Taurida"
replace region = "Tver'" if region == "Tver"
replace region = "Volynia" if region == "Volhynia"
replace region = "Yaroslavl'" if region == "Yaroslavl"

merge 1:1 region using ../work
gen serfdom = serfs/ruralpopn_1858
drop _m serfs
*Infer value for Ufa from value for Orenburg
list serfdom if region == "Orenburg"
replace serfdom = .0696855 if region == "Ufa"
save ../work, replace

*Soil quality
use ../work, clear
keep region
save ../temp, replace
use soils.dta, clear
rename name_gis region
replace region = "Nizhegorod" if region == "Nizhniy Novgorod"
merge m:1 region using ../temp
drop if _m ~= 3
drop _m russia1d_i
gen prcntsoilt = (soilareakm/areasqkm)*100
drop areasqkm soilareakm
collapse (sum) prcntsoilt, by(su_sym90 region)
replace su_sym90 = "OTHER" if su_sym90 == ""
reshape wide prcntsoilt, i(region) j(su_sym90, string)
mvencode prcntsoil*, mv(.=0) override
*Conmpare to FGO .do file: soil type VR not found in this sample
gen goodsoil = (prcntsoiltCH + prcntsoiltGR + prcntsoiltHS + ///
prcntsoiltKS + prcntsoiltPH)/100
gen blacksoil = prcntsoiltCH/100
drop prcnt*
merge 1:1 region using ../work
drop _m

*Provinces for which PG data were updated by Malyavskii
gen update_Mal = 0
replace update_Mal = 1 if region == "Astrakhan'" | region == "Saratov" | ///
region == "Bessarabia" | region == "Vitebsk" | region == "Vladimir" | ///
region == "Volynia" | region == "Voronezh" | region == "Vyatka" | ///
region == "Kazan'" | region == "Kiev" | region == "Kostroma" | ///
region == "Yaroslavl'" | region == "Kursk" | region == "Mogilev" | ///
region == "Moscow" | region == "Nizhegorod" | region == "Novgorod" | ///
region == "Orel" | region == "Penza" | region == "Perm'" | ///
region == "St.Petersburg" | region == "Pskov" | region == "Estland" | ///
region == "Podolsk" | region == "Ryazan'" | region == "Samara" | ///
region == "Simbirsk" | region == "Tambov" | region == "Tver'" | ///
region == "Tula" | region == "Ufa" | region == "Khar'kov"

*German occupation in 1917
gen occupied = 0
replace occupied = 1 if region == "Grodno" | region == "Kovno" | ///
region == "Kurland" |  region == "Minsk" | region == "Vilna" | ///
region == "Volynia" | region == "Estland" | region == "Lifland"

*Data labels for scatterplots
merge 1:1 region using region_short.dta
drop _m

*Changing St.Petersburg to Petrograd and other province name changes
gen region_17 = region
replace region_17 = "Petrograd" if region_17=="St.Petersburg"
replace region_17 = "Bessarabiia" if region_17 == "Bessarabia"
replace region_17 = "Don Host" if region_17 == "Don"
replace region_17 = "Nizhnii Novgorod" if region_17 == "Nizhegorod"
replace region_17 = "Podoliia" if region_17 == "Podolsk"
replace region_17 = "Riazan'" if region_17 == "Ryazan'"
replace region_17 = "Vil'na" if region_17 == "Vilna"
replace region_17 = "Volyn'" if region_17 == "Volynia"
replace region_17 = "Viatka" if region_17 == "Vyatka"
replace region_17 = "Iaroslavl'" if region_17 == "Yaroslavl'"
replace region_17 = "Ekaterinoslav" if region_17 == "Yekaterinoslav"
save ../work, replace


************************************
*****SPATIAL WEIGHTING MATRICES*****
************************************

cd ..

*Inverse-distance matrix
clear mata
shp2dta using GIS/SR_sample.shp, database(province) coordinates(provincexy) ///
genid(id) gencentroid(c) replace
use province, clear
rename NAME region
merge 1:1 region using work
drop _m
save work, replace
spmat idistance distmat x_c y_c, id(id)

*Contiguity matrix (produced by hand as the spmat contiguity 
*procedure produced results with some errors) 
spmat import contigmat using GIS/contigmat.txt  


*********************
*****DESCRIPTION*****
*********************

*Sumstats for Table 1 
sum PG_total_1917_pc MAL_totaltoOct_1917_pc  ///
PG_seizure_estate_as1917_pc PG_seizure_rest_as1917_pc ///
PG_comp_land_rent_as1917_pc PG_prevent_as1917_pc PG_other_as1917_pc ///
goodsoil serfdom occupied update 

*Correlations between PG and Malyavskii data
corr PG_total_1917_pc MAL_totaltoSept_1917_pc MAL_totaltoOct_1917_pc
corr PG_total_1917_pc MAL_totaltoOct_1917_pc
corr PG_total_1917_pc MAL_totaltoOct_1917_pc if update == 1

*Compare number of disturbances in updated and unupdated provinces
sum PG_total_1917_pc if update == 0
sum PG_total_1917_pc if update == 1


******************
*****ANALYSIS*****
******************

**
**Table 2
**

local vector "goodsoil serfdom occupied"

spreg gs2sls PG_total_1917_pc `vector', id(id) elmat(contigmat) het
outreg2 using table_SR_total, excel bdec(3) ctitle(PG) replace

spreg gs2sls MAL_totaltoOct_1917_pc `vector', id(id) elmat(contigmat) het
outreg2 using table_SR_total, excel bdec(3) ctitle(MAL)

spreg gs2sls MAL_totaltoOct_1917_pc `vector' update, id(id) elmat(contigmat) het
outreg2 using table_SR_total, excel bdec(3) ctitle(MAL)

spreg gs2sls PG_total_1917_pc `vector', id(id) elmat(distmat) het
outreg2 using table_SR_total, excel bdec(3) ctitle(PG)

spreg gs2sls MAL_totaltoOct_1917_pc `vector', id(id) elmat(distmat) het
outreg2 using table_SR_total, excel bdec(3) ctitle(MAL)

spreg gs2sls MAL_totaltoOct_1917_pc `vector' update, id(id) elmat(distmat) het
outreg2 using table_SR_total, excel bdec(3) ctitle(MAL)

**
**Regressions for Figure 4. Categories of disturbances. 
**PG data: April-September totals.
**
 
*Standardized variables
zscore PG_seizure_estate_as1917_pc PG_seizure_rest_as1917_pc ///
PG_comp_land_rent_as1917_pc PG_prevent_as1917_pc PG_other_as1917_pc ///
goodsoil serfdom occupied

local vector "z_goodsoil z_serfdom z_occupied"

*Estate seizures 
spreg gs2sls z_PG_seizure_estate_as1917_pc `vector', id(id) noconstant ///
elmat(contigmat) het
eststo estate_1
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(estate) replace
spreg gs2sls z_PG_seizure_estate_as1917_pc `vector', noconstant id(id) ///
elmat(distmat) het
eststo estate_2
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(estate)

*Other seizures
spreg gs2sls z_PG_seizure_rest_as1917_pc `vector', id(id) noconstant ///
elmat(contigmat) het
eststo seizure_1
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(seizure_other) 
spreg gs2sls z_PG_seizure_rest_as1917_pc `vector', id(id) noconstant ///
elmat(distmat) het
eststo seizure_2
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(seizure_other)

*Compulsory land rent
spreg gs2sls z_PG_comp_land_rent_as1917_pc `vector', id(id) noconstant ///
elmat(contigmat) het
eststo rent_1
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(rent)
spreg gs2sls z_PG_comp_land_rent_as1917_pc `vector', id(id) noconstant ///
elmat(distmat) het
eststo rent_2
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(rent)

*Prohibition of forest cutting and removal from work
spreg gs2sls z_PG_prevent_as1917_pc `vector', id(id) noconstant ///
elmat(contigmat) het
eststo prevent_1
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(prevent) 
spreg gs2sls z_PG_prevent_as1917_pc `vector', id(id) noconstant ///
elmat(distmat) het
eststo prevent_2
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(prevent)

*Other
spreg gs2sls z_PG_other_as1917_pc `vector', id(id) noconstant ///
elmat(contigmat) het
eststo other_1
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(other) 
spreg gs2sls z_PG_other_as1917_pc `vector', id(id) noconstant elmat(distmat) het
eststo other_2
outreg2 using table_SR_cat_AS, excel bdec(3) ctitle(other)


*****************
*****FIGURES*****
*****************

*Color scheme
set scheme s2mono

**
**Figure 3. Scatterplots of unrest, serfdom and fertile soil
**

*Changing position of province labels
gen pos = 3
replace pos = 10 if region_17 == "Petrograd"

gen pos_1 = 3
replace pos_1 = 9 if region_17 == "Petrograd" | region_17 == "Ekaterinoslav" ///
| region_17 == "Podoliia" | region_17 == "Kherson" | ///
region_17 == "Bessarabiia" | region_17 == "Khar'kov"
replace pos_1 = 6 if region_17 == "Voronezh"
replace pos_1 = 12 if region_17 == "Samara" | region_17 == "Kursk" ///
| region_17 == "Poltava" | region_17 == "Orenburg" | region_17 == "Petrograd"

graph drop _all

*Figure 3. Upper left panel.
graph twoway (lfit PG_total_1917_pc goodsoil, lwidth(0.7)) (scatter ///
PG_total_1917_pc goodsoil, mlabv(pos_1) mlabsize(small) ///
msymbol(O) mlabel(region_17)),  ///
xtitle({bf:Fertile soil}) ytitle({bf:Disturbances per capita}) ///
title({bf:Provisional Government data}, size(medlarge)) legend(off) ///
xscale(r(0.9)) xlabel(0 0.2 0.4 0.6 0.8) name(pgsoil) 
*graph export corr_PG_soil.png, width(7000) height(6000) replace   

*Figure 3. Upper right panel. 
graph twoway (lfit MAL_totaltoOct_1917_pc goodsoil, lwidth(0.7)) ///
(scatter MAL_totaltoOct_1917_pc goodsoil, mlabv(pos_1) mlabsize(small) ///
msymbol(O) mlabel(region_17)),  ///
xtitle({bf:Fertile soil}) ytitle(" ") title({bf:Malyavskii data}, ///
size(medlarge)) legend(off) xscale(r(0.9)) xlabel(0 0.2 0.4 0.6 0.8) ///
name(malsoil)
*graph export corr_MAL_soil.png, width(7000) height(6000) replace

*Figure 3. Lower left panel.
graph twoway (lfit PG_total_1917_pc serfdom, lwidth(0.7)) ///
(scatter PG_total_1917_pc serfdom, mlabv(pos) mlabsize(small) msymbol(O) ///
mlabel(region_17)), xtitle({bf:Serfdom}) ///
ytitle({bf:Disturbances per capita}) title(" ", size(medsmall)) ///
xscale(r(0.85)) xlabel(0 0.2 0.4 0.6 0.8) legend(off) name(pgserf)
*graph export corr_PG_serf.png, width(7000) height(6000) replace

*Figure 3. Lower right panel. 
graph twoway (lfit MAL_totaltoOct_1917_pc serfdom, lwidth(0.7)) ///
(scatter MAL_totaltoOct_1917_pc serfdom, mlabv(pos) mlabsize(small) ///
msymbol(O) mlabel(region_17)), xtitle({bf:Serfdom}) ytitle(" ") ///
title(" ", size(medsmall)) xscale(r(0.85)) xlabel(0 0.2 0.4 0.6 0.8) ///
legend(off) name(malserf) 
*graph export corr_MAL_serf.png, width(7000) height(6000) replace

*Combine the four graphs
graph combine pgsoil malsoil pgserf  malserf, col(2) imargin(vsmall) iscale(*.7)
*graph export correlations.tif, width(7000) replace
*graph export correlations.pdf, replace

**
**Figure 4. Coefplot with standardized coefficients
**

coefplot estate_1, bylabel(Model 1) || estate_2, bylabel(Model 2) || ///
seizure_1, bylabel(Model 1) || seizure_2, bylabel(Model 2) || ///
rent_1, bylabel(Model 1) || rent_2, bylabel(Model 2) || ///
prevent_1, bylabel(Model 1) || prevent_2, bylabel(Model 2) || other_1, ///
bylabel(Model 1) || other_2, bylabel(Model 2) ///
drop(_cons z_occupied) xline(0) bycoefs byopts(xrescale) mlcolor(none) /// 
ciopts(lwidth(*2.5)) msize(2.5) levels(95) baselevels  ///
headings(1 (2) = "{bf:Estate} {bf:seizures}" ///
3 (4) = "{bf:Other} {bf:seizures}" ///
5 (6) = "{bf:Compulsory} {bf:land} {bf:rental}" ///
7 (8) = "{bf:Disruption} {bf:of} {bf:activities}" 9 (10) = "{bf:Other}") ///
coeflabels(z_goodsoil = "{bf:Fertile soil}" z_serfdom="{bf:Serfdom}")
*graph export coefplot.png, width(7000) height(6000) replace
*graph export coefplot.tif, width(7000) height(6000) replace
*graph export coefplot.pdf, replace
